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Abstract 

We consider the solid-solid interactions in the two body 
problem. The relative equilibria have been previously 
studied analytically and general motions were numerically 
analyzed using some expansion of the gravitational poten- 
tial up to the second order, but only when there are no 
direct interactions between the orientation of the bodies. 
Here we expand the potential up to the fourth order and 
we show that the secular problem obtained after averag- 
ing over fast angles, as for the precession model of Boue 
and Laskar [Boue, G., Laskar, J., 2006. Icarus 185, 312- 
330] , is integrablc, but not trivially. We describe the 
general features of the motions and we provide explicit 
analytical approximations for the solutions. We demon- 
strate that the general solution of the secular system can 
be decomposed as a uniform precession around the to- 
tal angular momentum and a periodic symmetric orbit 
in the precessing frame. More generally, we show that 
for a general n-body system of rigid bodies in gravita- 
tional interaction, the regular quasiperiodic solutions can 
be decomposed into a uniform precession around the total 
angular momentum, and a quasiperiodic motion with one 
frequency less in the precessing frame. 
Keywords : CELESTIAL MECHANICS, SOLID-SOLID 
INTERACTION, BINARY ASTEROID DYNAMICS, 
ROTATIONAL DYNAMICS 



1 Introduction 

We consider here two rigid bodies orbiting each other. 
The main purpose of this work is to determine the long 
term evolution of their spin orientation and to a lower 
extent, the orientation of the orbital plane. Examples 
of such systems are binary asteroids or a planet with a 
massive satellite. 
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If the two bodies are spherical, then the translational 
and the rotational motions are independent (e.g. Du- 
boshin, 1958). In that case, the orbit is purely keplerian 
and the proper rotation of the bodies arc uniform. Gen- 
eral problems with triaxial bodies are more complicated, 
and usually non integrable. Even formal expansions of the 
gravitational potential or the proof of their convergence 
can be an issue (Borderies, 1978; Paul, 1988; Tricarico, 
2008). In some cases, especially for slow rotations close to 
low order spin-orbit resonances, the spin evolution of rigid 
bodies of irregular shape can be strongly chaotic (Wisdom 
et al., 1984; Wisdom, 1987), but we will not consider this 
situation in the present paper where we focus on regular 
and quasiperiodic motions. 

Stationary solutions of spin evolution are known in 
the case of a triaxial satellite orbiting a central spheri- 
cal planet (Abul'naga and Barkin, 1979). In their paper, 
Abul'naga and Barkin used canonical coordinates, based 
on the Eulcr angles, to set the orientation of the satellite. 
On the contrary, in 1991, Wang et al. also studied rela- 
tive equilibria but with a vectorial approach that enabled 
them to analyze easily the stability of those solutions. For 
a review of different formalisms that can be used in rigid 
body problems, see (Borisov and Mamacv, 2005). 

The vectorial approach turned out to be also powerful 
for the study of relative equilibria of two triaxial bodies 
orbiting each other (Maciejewski, 1995). General motions 
of this problem were studied by Fahnestock and Scheeres 
in 2008 (hereafter FS08) in the case of the typical binary 
asteroid system called 1999 KW4. For that, the authors 
expanded the gravitational potential up to the second or- 
der only. In this approximation, there is no direct interac- 
tion between the orientation of the two bodies. Ashenberg 
gave in 2007 the expression of the gravitational potential 
expanded up to the fourth order but didn't study the so- 
lutions. 

In (Boue and Laskar, 2006) (hereafter BL06) we gave a 
new method to study the long term evolution of solid body 
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orientations in the case of a star-planet-satellite problem 
where only the planet is assumed to be rigid. This method 
used a similar vectorial approach as Wang et al. (1991), 
plus some averaging over the fast angles. We showed that 
the secular evolution of this system is integrable and pro- 
vided the general solution. 

In the present paper, we show that the problem of two 
triaxial bodies orbiting each other is very similar to the 
star-planet-satellite problem and thus can be treated in 
the same way. 

In the section 2, we compute the Hamiltonian govern- 
ing the evolution of two interacting rigid bodies. The 
gravitational potential is expanded up to the fourth or- 
der and averaged over fast angles. The resulting secular 
Hamiltonian is a function of three vectors only: the or- 
bital angular momentum and the angular momenta of the 
two bodies. 

In a next step (section 3), we show that the secular 
problem is integrable but not trivially (i.e. it cannot be 
reduced to a scalar first order differential equation that 
can be integrated by quadrature). The general solution 
is the product of a uniform rotation of the three vectors 
(global precession around the total angular momentum) 
by a periodic motion (nutation). We prove also that in a 
frame rotating with the precession frequency, the nutation 
loops described by the three vectors are all symmetric 
with respect to a same plane containing the total angular 
momentum. We then derive analytical approximations 
of the two frequencies of the secular problem with their 
amplitudes. These formulas need averaged quantities that 
can be computed recursively. However we found that the 
first iteration already gives satisfactory results. 

In section [5j we consider the general case of a n-body 
system of rigid bodies in gravitational interaction, and 
we demonstrate that the regular quasiperiodic solutions 
of these systems can, in a similar way, be decomposed into 
a uniform precession, and a quasiperiodic motion in the 
precessing frame. 

Finally, we compare our results with those of FS08 
on the typical binary asteroid system 1999 KW4. We 
show that their analytical expression of the precession fre- 
quency corresponds to the simple case of a point mass or- 
biting an oblate body treated in BL06. We then integrate 
numerically from the full Hamiltonian, an example of a 
doubly asynchronous system where the FS08 expression 
of the precession frequency does not apply. We compare 
the results with the output of the averaged Hamiltonian 
and with our numerical approximation and show that they 
are in good agreement. 



2 Fundamental equations 

We are considering a two rigid body problem in which 
the interaction is purely gravitational with no dissipa- 
tive effects. Let m\ and m 2 be the masses of the two 
solids. Hereafter the mass m 2 is called the satellite or the 
secondary and the mass ui\ the primary. It should be 
stressed that this notation does not imply any constraint 
on the ratio of the masses which can even be equal to one. 

The configuration of the system is described by the 
position vector r of the satellite barycenter relative to 
the primary barycenter and their orientation expressed in 
an invariant reference frame. The orientations are given 
by the coordinates of the principal axes (Ji, Ji, K{) and 
(i 2, J 2, K 2 ) in which the two inertia tensors, respectively 
Ii and I 2 , of the primary and of the secondary are diag- 
onal [Ji = diag(Ai, Bi, C\) and 1 2 = diag(A 2 , B 2 , C 2 )]. 

The Hamiltonian of this problem can be split into 



h = h 1 



Hp 



Hr 



(1) 



where Ht is the Hamiltonian of the free translation of the 
reduced point mass (3 = m 1 m2/(mi + m 2 ), He describes 
the free rigid rotation of the two bodies and Hi contains 
the gravitational interaction. 

The Hamiltonian of the free point mass is 



Ht — 



2/3 



(2) 



where f = /3r is the conjugate momentum of r. 

Let Gi and G 2 be respectively the angular momentum 
of the primary and of the satellite. The Hamiltonian of 
the free rotation is 



H h 



t G\T 1 1 G\ 



t G 2 T 2 1 G 2 



(3) 



where the superscript * in *x or 1 A denotes the transpose 
of any vector x or matrix A. It can be expressed in terms 
of the principal bases of the two bodies as follows 



H E = 



(GW1) 2 , (G 1 -J 1 ) 2 , (G 1 -K 1 f 



2A 1 2B 1 2d 

, (G 2 -I 2 f , (G 2 -J 2 ) 2 , (G 2 -K 2 f 



2A 2 



2B 2 



2C 2 



(4) 



The interaction between the two solid bodies is the fol- 
lowing double integral 



H, 



Q dmi dm 2 
|r + r 2 -n| 



(5) 



where ri and r 2 are respectively computed relative to the 
primary and satellite barycenters (cf Fig. [IJ and describe 
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Figure 1: Coordinates definition. 



the two volumes. This part of the Hamiltonian can be 
expanded in terms of Legendre polynomials and will be 
written as a function of (r, Ji, Jx, K\, I2, J2; K 2 ) in the 
section 12.31 



2.1 Equations of motion 

The full Hamiltonian is written in the non-canonical co- 
ordinates (r, f , 1 1, J 1, Ki, Gi, J 2j J 2, K 2 , G 2 ). Thus, al- 
though the components (r,f) keep the standard sym- 
plectic structure, (Jx, «7i, K x, Gi) on the one hand and 
(J2, J 2, K 2 , G2) on the other hand possess the Euler- 
Poisson structure which leads to the following equations 
of motion (Borisov and Mamaev, 2005) 



2.3 Gravitational potential 



The distance between the two bodies is assumed to be 
large in comparison to their size. Thus, in the expres- 
sion of the gravitational potential p\ = ||ri||/||r| 
and pi — ||r2|j/||r|| are two small parameters. It can 
then be expanded in terms of Legendre polynomials (see 
Appendix A). As it is shown below (equation 13), the 



expansion up to the second order does not contain any 
interaction due to the relative orientation of the bodies. 
We thus choose to expand the gravitational potential up 
to the fourth order. In the computation appear integrals 
such as J drrii or J r^Yidrrii, i = 1,2 which can be 
expressed in terms of moments of inertia 



r = V~ r H, r = -V r H 

G = V/W x I + VjH x J + v K n x K + V G H x G 
i = V G H x 7, j = V G H x J, K = V G H x K. 

(6) 

We choose these non-canonical coordinates instead of 
symplectic ones because of the simplicity of the resulting 
equations which already resemble equations of precession. 



2.2 First simplification 

In the previous paragraphs, the Hamiltonian contains the 
three vectors of the principal frame (J, J, K) of each 
body. Nevertheless, only two vectors per solid are nec- 
essary insofar as the third can be expressed as the wedge 
product of the other two. We choose to keep J and K. 

The Hamiltonian of the free rotation of the two rigid 



J r i dm i = 2 ' 

Jr^ridm,^ — — + — Id 

+ (B, - A l )I l t I l + (Bi - C^K/Ki , 

(8) 

with Id beeing the identity matrix in M 3 . But higher de- 
gree integrals such as J rf drrii also appear. To compute 
these integrals, one needs more information about the 
bodies. However, moments of inertia are already hardly 
known, at least for satellites. It is thus not relevant to 
add new unconstrained parameters. But such integrals of 
inertia can be expressed as functions of Ai, Bi, Ci assum- 
ing that the bodies are homogeneous ellipsoids. Indeed, 
let (xi, yi, Zi) be the coordinates in the principal frame of 
a running point of the body i, and Ip, q ,r-i — J x^y^zl drrii 
be its integrals of inertia. Because of the three symmetry 
planes of homogeneous ellipsoids, I p , q , r :i vanishes when- 
ever one of p, q, r is odd. Thus all the third order integrals 
of inertia cancel, and the only non zero fourth order inte- 
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grals of inertia are (see Appendix B) 



order terms expression is classical and given by 



xf dm; 



yf drrii = 
zf drrii = 
y?zfdmi = 
z 2 x 2 drrii = 
x 2 y 2 drrii = 



15 

28m., 
15 

28m., 
15 

28m, 
5 

28m., 
5 

28m., 
5 

28m., 



(-Ai + B t + Ci) 2 ■ 

{Ai -Bi + d) 2 ■ 

(A % + Bi — Ci) 2 ■ 

{Ai -Bi + C i )(A l + Bi- d) ■ 

{Ai + Bi-dX-Ai + Bi + d) 

{-Ai + Bi + CJiAi-Bi + Ci) 



(9) 

In search of generality, we now forget the assumption 
of homogeneous ellipsoids. We only keep the symmetry 
plane hypothesis that cancels odd integrals. Setting 



H i 2) =~\% t m i(^2 - 2B 2 + C 2 ) + m 2 (A 1 - 2B X + C x )] 
2 r J 



3 Qm\ 



2 r 3 

3 Qm 2 
2^ 



(B 2 - A 2 ) (u- J 2 ) 2 + (B 2 - C 2 ) (u- K 2 f 



(B 1 - A,) (u- Ji) 2 + (B 1 - C x ) (u- K x f 



(13) 



X >=J 


[ x\ drrii 


p > = ) 






«-J 


yf drrii 




f z 2 x 2 drrii 5 


(10) 


Z '=J 


zf drrii 


R ' = J 


f x 2 y 2 drrii 





where u = r/r is the direction vector of r. As mentioned 
before, this expression does not contain body-body in- 
teractions but only spin-orbit ones such as (u- Ki) 2 or 
(u • K 2 ) 2 . The fourth order terms expression is given 



the integrals appearing in the expansion of the gravita- 
tional potential become 



r'i dm; 



J (s- i-j) 4 drrii 



in (14 1. In contrast to the second order terms, among 
the fourth order terms there are direct interactions be- 
tween the two orientations such as (-KV K 2 ) 2 . A similar 
expression was published recently in (Ashenberg, 2007). 
Although more terms are present in Ashenberg's paper 
because we have made here the additional assumption of 
symmetry of the rigid bodies, we could compare our ex- 
pression successfully with the one of Ashenberg, except 
for a difference in a coefficient that may come from a mis- 
print in Ashenberg's papeiQ 
2s 2 [(3R l - Yi)(s- 1,) 2 + (3P, - Yi)(s- K t ) 2 ] 

2[Yi - 3(P< -Q, + Ri)](s- Ii) 2 {s- K,) 2 ; 



X. t + Y l + Z l + 2P, + 2Q l + 2Ri 

YiS 4 + (Xi + Yi— 6Ri)(s- Ji) 4 
+{Z l + Y l -QP i ){s- Ki) A 



rfvi'Ti drrii 



(Yi + Ri + Pi)Id + {Xi -Yi + Qt- P t )I l t I l 
+(Zi -Y + Q t ~ R i )K l t K l 

(11) 

where s is any vector and i = 1,2. 

With these results, the expansion of the potential gives 
the zeroth order term 



hP=-^- (12) 

v 

1 In (Ashenberg, 2007), there is a misprint in the ex pre ssion of 
Vgjg/, Eq. (20). The coefficient -35/(4r 5 ) in Eq. Jill of the 
where (J, = Q(mi + m 2 ). This is the well known gravita- curren t paper corresponds to a coefficient -G/(SR 5 ) in Ashenberg's 
tional interaction between two point masses. The Second notations whereas it is written -G/(5R 5 ) in (Ashenberg, 2007). 
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3 Q 



I 4 r 5 



(A 2 - 2B 2 + C 2 )(A 1 - 2Bi + Ci) 



-ma[Xi + -Yt +Z-L- 8P1 + 2Q 1 - 8P1] 

l - mi [X 2 + ^Y 2 + Z 2 - 8P 2 + 2Q 2 - 8P 2 ] 

2(B 1 -A 1 )(B 2 -A 2 )(I 1 -I 2 ) 2 
2(B 1 -A 1 )(B 2 -C 2 )(I 1 -K 2 ) 2 
2{B 1 -C 1 ){B 2 -A 2 ){K 1 -I 2 ) 2 
2{B 1 ~C l ){B 2 -C 2 ){K 1 -K 2 f 

5(A 2 -2B 2 + C 2 )(B 1 - Ax) 
20 

- to2(5X! + —Y 1 - 5P X + 5Qi - 35jRi) (u- Iif 

b{A 2 -2B 2 + C 2 ){B 1 ~C 1 ) 
20 

- m 2 {hZ l + — Y l - 35Pi + 5Qi - 5Pi) (u- 

5(A 1 -2 J B 1 + Ci)(B 2 -^ 2 ) 
20 

- mi(5X 2 + y7 2 - 5P 2 + 5Q 2 - 35i? 2 ) (u- I 2 y 

5(A 1 -2B 1 + C 1 ){B 2 -C 2 ) 
20 

- mi (5Z 2 + —Y 2 - 35P 2 + 5Q 2 - 5i? 2 ) (u- K 2 y 



-20 

35 
T 

35 
¥ 

-35 



(Pi - Ai)(u- + (Pi - Ci)(u- -ftTi)Ki 
(P 2 - A 2 )(u- I a )I 2 + (P 2 - C 2 )(u- K 2 )K 2 
m 2 \{X 1 + Y x - 6Pi ) (u- 1 1 ) 4 + (Zi + Yi - 6Pi ) (u- X x ) 4 
+ 2(Fi - 3Pi + 3Qi - 3Pi)(u- /i) 2 (u- fd) 21 



-mi 



(x 2 +y 2 -6P 2 )(u-/ 2 ) 4 +(z 2 +r 2 -6P 2 )(u-x 2 ) 4 

2(F 2 - 3P 2 + 3Q 2 - 3P 2 )(u- J 2 ) 2 (u- ii: 2 ) 21 



(Pi - At)(u- Ji) 2 + (Pi - Ci)(u- Ki) 2 
(P 2 - A 2 )(u- I 2 f + (P 2 - C 2 )(u- 1<: 2 ) 2 



(14) 



The full Hamiltonian (ff) [2| [l2j [l3| an d [l4|> together with 
the equations of motion (cf section 2.1 1 enable the integra- 
tion of the system. The evolution of this system contains 
fast motions like the rotation of each body around their 
axis or the orbital revolution. In comparison, the two spin 
axes as well as the orientation of the orbital plane undergo 
secular evolutions. In the following, fast motions are av- 




Figure 2: Definition of Andoyer's coordinates. (i,j,k) is a fixed 
reference frame, and (I, J, K) the reference frame of the principal 
axis of inertia of the solid body. The Andoyer action variables are 
(G, H = G ■ k, L = G ■ K) with the associated angles (g, h, I) (An- 
doyer, 1923). 



eraged in the purpose of studying the long term evolution 
only. 

2.4 Averaging 

In this section, we average the Hamiltonian independently 
over all fast angles: proper rotations and orbital motion. 
Although this method is strictly valid for non resonant 
cases only, we will show (in section [6]) an application to a 
typical primary-asynchronous, secondary-synchronous bi- 
nary asteroid system where the motion is regular. The 
method still gives very acceptable results. In the follow- 
ing, we forget the subscripts 1 and 2 whenever we consider 
any of the two bodies without distinction. 

To average over proper rotations, Andoyer variables 
(G,H,L,g,h,l) as described in Fig. ([2| are well suited. 
In a first step, the dependency of the full Hamiltonian on 
1 1 and I 2 is removed by averaging over l\ and l 2 . We 
have 

( cos I 

(15) 

(n,ti',K) 

where n is defined in Fig. Q and n' = K x n. The 




n! and K 


are independent of I, thus 


Ml 


= 0; 




= l{Id-K*K) ; 


<(-') 4 >, 


= 3 -[s 2 -(s-K) 2 ] 2 , 



(16) 
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where s is again any vector. After this averaging, the 
Hamiltonian of the free rotation becomes 



2A'x + I d A' i 



2A'- 



1_ J_\{G 2 -K 2 f 

C 2 A' : 



where 



11/1 1 

A* ~ 2\A + B 



(17) 



(18) 



And the second and the fourth order terms of the inter- 
action 



H 



(2) 



QCim 2 
2r 3 



QC 2 rn 



2r 3 



' [l-3(u-K 2 ) 2 ], (19) 



vectors Ni, Ui and w are independent of g, so 



= (cos J)w 



(K t K) g = i( s in 2 J)7rf+ M - ^ sin 2 j) w'w 



{(s-K) 4 ) = (l-5sin 2 J+^sin 4 j) (s-w) 



-3 sin 2 J ( 1 - | sin 2 J ) s 2 (s- w) 2 



3 „:„4 



• 4 t 4 

sin Js , 
8 

(23) 

where s is any vector. After averaging over g\ and <? 2 , the 
conjugated momenta G\ and G 2 become constant. The 
averaged Euler Hamiltonian which depends only on G\ 
and G 2 



H 



3 Qm 2 T>i 



3 Gm 1 T> 2 
8 

3 ^CiC 2 



l-10(u-K 1 ) 2 + y(u-K 1 ) 4 
l-10(u-K 2 ) 2 + y(u-X 2 ) 4 



l + 2(fsT 1 -i<: 2 ) 2 -5(u- J Fs: 1 ) 5 



4 r 5 

-5(u- K 2 f - 20(u- lfi)(u- K^K^K^ 
+35(u-K 1 ) 2 (u-K 2 ) 21 



where C = C - ( A + B)/2 and 



(20) 



V=l(X + Y) + Z -3(P + Q)+ 3 I R. (21) 
8 4 



In a next step, the averaging over the angle g is per- 
formed. This corresponds to the averaging of K around 
w = G/G (cf Fig. [2]). Indeed, in the general case the an- 
gular momentum G is not aligned with the axis of max- 
imum inertia K, which is implicitly assumed in the gy- 
roscopic approximation. Instead, if there is an angle J 
between these two vectors then 



sin J sin g 
K = I — sin J cos g 

cos J 



(22) 



(Ni,Uj.,w) 

where Ni is defined in Fig. [2] and Ui = w x Ni. The 



(He) 



f cos 2 Ji sin 2 Ji \ G\ 



l,g 



A'x 



cos 2 J 2 sin 2 J 2 \ G 2 
C 2 + A' 2 ~ 



(24) 



is now a constant and can be ignored. In this expression, 



A' is still the harmonic mean of A and B (18 1. The only 



change in the interaction is the substitution of C and V 
in ppt by 



C" = ( 1 - I sin 2 J j I 7 



£>' = 1-5 sin 2 J 



35 



sin 4 J V 



(25) 



and (Ki,K 2 ) by (wi,w 2 ). For fast rotating non-rigid 
bodies, the angle J is assumed to be very small as a result 
of internal dissipation (J ss 1CP 7 radians for the Earth). 
In that case, the gyroscopic approximation J = is a 
good approximation since the correction obtained after 
averaging over fast angles is in 0(sin 2 J). Nevertheless, 
for slow rotating triaxial asteroids, the angle J may be 
large and the gyroscopic approximation may not be valid. 

In a third step the Hamiltonian is averaged over the 
orbital motion. First over the mean anomaly M, and 
then over the longitude of periapse uj. The first average is 
computed using the formulas of the Appendix C and for 



(i 



the second one, we have similar equations as ( 16 1 



where 



k 3 m 2 C[ ~ -h{C[C 2 + m 2 V 1 ) 



; 



(26) 



<(-') 4 >. 



'[* 2 -(-w) 2 ] 2 



where I now denotes direction of the periapse and w the 
normal of the orbit. The resulting secular Hamiltonian 



H s ~ VHl u l 2 , gi ,g 2 ,M,u 1S 



thus 



( cos 2 Ji sin 2 J\ \ G\ ( cos 2 J 2 sin 2 J 2 \ G 



2a 



Co 



A' 



[m 2 C[(l - 3x 2 ) + mi C' 2 {\ - 3y 2 )] 
3 



4a 3 (l-e 2 ) 3 / 2 

9 G ( 

32a 5 (l-e 2 ) 7 /2 { + 2 6 

C{C' 2 {\ - bx 2 - by 2 + 2z 2 - 2Qxyz + 35x 2 y 2 ) 
+ ^fl (l - 10, 2 + fx 4 

+^ (l-l0y 2 + 3 p 

(27) 

where x = (w- wj, y = (w- W2) and z = (wi • W2). Let 
us write H s in the more compact form 



tt a 2 b 2 

H, = --X —-V - 



-z +dxyz— —x 



La 



x z y z + f) 
(28) 



/c 3 miC 2 - -k 5 (C[C' 2 + niiV' 2 ) 



k$CiC 2 



= 



35 

¥ 

35 
TT 

35 



-fc 3 (m 2 C; + miC 2 ) - -fc 5 (2C;C 2 + m 2 D; + mi ^) 



with 



(He) 

h 

h 



l <9 2a 
3 Q 



2a 3 (l-e 2 ) 3 / 2 

9 G 

8a 5 (l-e 2 ) 7 /2 



1 + -e 4 



(29) 
(30) 



3 Secular equations 



The secular Hamiltonian H s (28) is similar to the one 



obtained in BL06 although its expression is slightly more 
complicated. The difference with BL06 is that the secular 
Hamiltonian is not anymore the equation of an ellipsoid in 
(x, y, z). A few results in BL06 were proved for this special 
surface. We recall here the main steps of the derivation 
of the solutions adapted to the new surface defined by the 
current secular Hamiltonian. 

The Hamiltonian H s is only a function of the angular 
momenta (G, Gi, G 2 ). The equations of motion of these 
quantities are 



G = V G H S x G , 
G 2 = V G . 2 H S x G 2 



(31) 



We thus have G G = d d 
that the norms 7 = ||G||, /3 - 



G 2 - G 2 = which means 
IIGill and a = ||G 2 || are 



7 



constant. It is thus possible to write the general equations 
directly in terms of the unit vectors (w, Wi, W2) 



y = w- w 2 and z = wi - w 2 . From the equations (33 1, one 
can derive the new equations of motion 



w = -V W H S x w , 

7 

wi = ^V Wl H s x wi 

P 

w 2 = —V vf2 H s x w 2 



(32) 



From the expression of the secular Hamiltonian ( 28 ) , we 
get 

P 1 
w = Wi x w w 2 x w , 

7 7 



where 



p s 

Wi = — — w x wx — -W2 X Wi , 

P P 

q s 

W 2 = W X W 2 Wi X W 2 , 

a a 



p = ax — dyz + cx 3 + gxy 2 , 
q = by - dxz + fy 3 + gx 2 y , 
s = cz — dxy . 



(33) 



The problem has 9 degrees of freedom, the coordinates 
of G, Gi and G 2 , and the equations ( 33p4 1 are non- 
linear. At first glance the resolution is difficult. There 
are 7 first integrals 

llwll = 1 



Wl 



1 



o.r 



Wall = 1 
by 2 - 
/3wi 



7W 



cz — 2Dxyz ■ 
- arw 2 = Wq 



2 y 



Qx 2 y 2 = -2H S 



(35) 

where Wq is the total angular momentum. Thus one 
misses one constant of motion to integrate the problem 
by quadrature. The next section shows how to solve the 
relative motion of the three vectors that contains enough 
constants of motion. 

3.1 Relative solution 

In the previous section, we have shown that the number of 
first integrals is not large enough to solve the full problem. 
But the number of degrees of freedom can be decreased by 
considering only the relative distance between the vectors. 
These distances are given by the dot products x = w-wi, 



V = 



g 
7 
s 
a 

P 

1 



(36) 



where v = (w x Wi) • w 2 is the volume defined by the 3 
vectors. It can be expressed in terms of x, y and z through 
the Gram determinant 



x 
1 

z 



z + 2xyz. 



(37) 



This restricted problem has only 3 degrees of freedom 
and 2 first integrals 



ax 



by 2 



cz 



2dxyz 

(34) 7/3a; + ajy + a/3z = K 



f 4 

2 y 



Qx 2 y 2 



-2H* 



(38) 

the second integral being simply derived from 2K = 
ll^ol| 2 — ( r y 2 +P 2 +a 2 )- The motion in (x, y, z) is thus inte- 
grable and the solution evolves in the intersection ^ of the 
quartic H s — Cte and the plane K — Cte [^] Moreover, 
the evolution is limited to the interior of the v 2 (x, y, z) = 
surface that will be henceforth called the Cassini berlin- 
got[^]as in BL06 (cf Fig. Outside this surface we would 
have v 2 < which is not possible (see BL06). 

3.1.1 Shape of the quartic surface 

The constraint H s = Cte defines a quartic surface Q in 
(x, y, z). Quartic surfaces can have very different shapes, 
nevertheless setting z' = 



2H, = ax 2 



^xy, one obtains 



" ' ~2 yA ~{ ~\) x2y2 

(39) 

which is a biquadratic. The new surface Q' defined by 



(39) is thus symmetric in x, y and z' . In (x 2 ,y 2 ,z') the 



surface Q! can be cither an ellipsoid, a paraboloid or an 
hyperboloid depending on the sign of 



ef 



9 

c 



(40) 



2 In the whole paper, Cte means any constant value. 

3 A berlingot is a famous tetrahedron hard candy with rounded 
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Figure 3: The surface v 2 (x,y, z) = 0. As v 2 > 0, the allowed space 
is the interior of this Cassini berlingot shaped volume. 



If 6 > then it is an ellipsoid and x, y, z' and thus z 
are bounded. In the other case, Q' is cither an elliptic 
paraboloid if 6 = or an hyperboloid of one or two sheets 
depending on the value of H s if 5 < 0. Thus, x, y, z are 
unbounded. 



CD 

Sh 

hC 

CD 





H 



From the definition of the coefficients a-Q ( 29 1 , S can 
be rewritten in the following form 
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J\ (degree) 

Figure 4: Shape of the surface Q' as a function of the angles J\ 
and J2 between the angular momentum and the axis of maximum 
inertia of each bodies. £ and Ti stand for ellipsoid and hyperboloid 
respectively. The solid curve in red delineates the £ and the Ti zones 
for r) = 25/9. The dashed blue curve corresponds to r\ = 25/4. See 
text for details. 





6 = 



y) klcfC?. (-11) 



Using the definition of the coefficients C and V ( 25 1 and 



(21 1, we get 

S = 



15 



h 2 c 2 r 2 

ft '5 L 'l l -'2 



7/1-5 sin 2 Ji 



35 



sin J\ 



X 1 - 5 shr J 2 



35 



sin J s 



(42) 



l-^sin 2 Jij n-^sin 2 J 2 



where r\ is a positive parameter related to the shapes of 
the rigid bodies 



V 



7\ mi2?im 2 2?2 

9) cici 



(43) 



Let us look to the range of the possible values of r\ in the 
case of an homogeneous ellipsoids. We have the relation 



between V and C given by the equations (|9| 



V 



15 
7m 



C 2 + \{B-Af 



(44) 



The lowest value of rj is thus obtained for A — B, i.e. for 
axisymmetric bodies. In that case, r/ m i n = 25/9. Con- 
versely, the largest value of 77 is attained when (B — A) 2 
is maximal, thus when B = C and A = 0, that is, in 
the limiting case where the bodies are extremly thin rods. 
In this second case, r) max = 25/4. So, for homogeneous 
ellipsoids, rj is constrained between ?7 m in and rj max . 

Figure [4] shows the domains where the surface Q! is an 
ellipsoid (£) or an hyperboloid (TL) as a function of the 
angles J\ and J 2 . The two sets of curve correspond to 
V = Vmin an d fj = ?7max- As 8 is a function of sin 2 Ji 
and sin 2 J 2 , the figure can be extended up to 180 degrees 
applying axial symmetry around the axis J± = 90 degrees 
and J 2 = 90 degrees. 

3.1.2 Description of the solutions 

In BL06, we show that when the surface Q is an ellip- 
soid then the evolution of (x, y, z) presents two kinds of 
solutions. We have called special solutions the solutions 
where is totally included in the Cassini berlingot B. 
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Figure 5: The shaded area corresponds to the region where v 2 > 
0, inside the Cassini berlingot B. The orbit intersects the Cassini 
bcrlingot B in t = t+ and t = r_ . 



This means that the vectors w, wi and W2 are never 
collinear. This happens only when the three vectors are 
almost orthogonal. The second class of solutions are the 
general solutions, more frequent in astronomical prob- 
lems, for which crosses the Cassini berlingot (Fig. |5|. 
In that case M = (x,y,z) does periodic returns inside 
the Cassini berlingot up to its surface and the volume v 
defined by the three vectors w, wi and w 2 is conserved 
over one period. In both cases, solutions are periodic. 

There are also special cases that happen when the orbit 
of (x, y, z) is tangent to the Cassini berlingot. At the 
tangency there is indeed a fixed point. In that state, the 
three vectors remain in a plane that precesses in time. It is 
called a Cassini state (Colombo, 1966; Peale, 1969; Ward, 
1975; BL06). If an initial condition is chosen along such 
special orbits but strictly outside the fixed point, then the 
system cannot reach the stationary point in finite time 
and it is the only case where x, y, z are not periodic. 

Here, we have the same results except when the quartic 
Q is unbounded. In that case, we cannot have special 
solutions. 

3.2 Global solution 

Knowing that x, y, z are periodic functions of time, it 
is possible to get general features of the global motion. 



For that, let us rewrite the secular equations (33 1 in a 



new form so as to obtain a linear differential system with 
periodic coefficients. 

Let us assume as in BL06 that the vectors (w,w 1 ,w 2 ) 
are not coplanar (v ^ 0). Let W be the ma- 
trix (w, Wi, W2) and V the Gram matrix of the basis 
(w,wi, w 2 ) 



V = 




(45) 



Using the expression of the wedge product in the basis 
(w, Wi,W2) (see the appendix B of BL06), the equations 



of motion ( 33 ) can be written in the following form 

W = vV^WA. (46) 

Here we correct a mistake^ in the demonstration of the 
proposition 1, given in section 4 in BL06 (see the erratum 
Boue and Laskar, 2008). 

In (46 1, vV~ x and A are matrices depending only on 



(x, y, z) that are periodic functions of period T. Indeed 




A 



_ <± 

7 

V 

7 







_ p 





and 



V~ 




(47) 



(48) 



Thus, if W(t) is a solution of ([46} , then W(t + T) is 
also a solution. Let us denote 



TZ T (t) = W(i + T)W(t) _1 



(49) 



We need to prove that TZr(t) is constant with t. As the 
Gram matrix V of the vectors (w(i), Wi(t), W2{t)) is T- 
periodic, the norm is conserved by linear transformation 
n T {t) that send W(t) into W(t + T), and K(t) is thus 
an isometry of K 3 . Moreover, this isometry is positive, as 
the volume v is conserved over a full period T (see section 
3.1.21. The invariance of the total angular momentum 



Wq (35 1 then implies that IZrit) is a rotation matrix of 
axis Wo. 

As IZrit) is a rotation in M 3 , we have for all Wj,Wj in 
{w,wi,w 2 }, 

w t (t + T) x Wj(t + T) = (n T {t)wi(t)) x (1Z T (t)™ 3 (t)) 
= K T (f)(w,(()xw J (i)) . 

(50) 



From the equations of motion ( 33 1 , we can thus derive 



W(t + T) =K T (t)W(t) 



(51) 



n T {t)W{t) (49), we 



(52) 



On the other hand, as W{t + T) 
deduce that for all t, 

n T {t)w{t) = 

TZrit) is thus a constant matrix TZx- Now, let us de- 
note TZ(t) the rotation of axis Wo and angle tOx/T (i.e. 
TZ(T) =TZ T ). We have 

4 In BL06, we have incorrectly stated that the averaged differen- 
tial system (51) could be written as VV = ~WB where B = vV —1 A is 
a matrix depending only on (x, y, z). In fact the correct expression 
is VV = vV -1 )^^-. In BL06, the proof following the equation (51) 
has to be modified. This is done in the present paper. The results 
remain identical. 
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Proposition 1 The complete solution VV(t) can be ex- 
pressed on the form 



This final expression is an explicit function of (ar, y, z). 
Indeed, from (35 1, one has 



W(t) = K(t)W(t) , 



(53) 



where W(t) is periodic with period T , and lZ(t) a uniform 
rotation of axis Wo and angle tOx/T. The motion has two 
periods: the (usually) short period T and the precession 
period 

T' = ^T. (54) 

3.3 Properties of the solution 

A more precise result on the periodic loops can be proved. 
But before, one needs to write the instantaneous preces- 
sion speed as a function of (x, y, z). 

3.3.1 Instantaneous precession rate 

Let us write the time derivative of the precession angle 
of w as a function of (x,y,z). The expressions for the 
other vectors can be obtained in the same way. The 
following approach is highly inspired by BL06. We set 
Wo = 1 1 Wo 1 1 the norm of the total angular momentum 
and w = Wq/Wq its direction vector. With 

C = w-w , (55) 

the projection L of w on the plane orthogonal to w is 

L = w C w o- (56) 

Assuming w ^ w , we get £ < 1. With L = ||L||, the 
expression of L gives 



and 



L = -- 



cc 



Moreover, setting £ = h/L, we get 

L = LI and t = L£ + 9(w Q x L) 
which yields to 

L 2 ~ L 2 + 9 2 (w x L) 2 — L 2 + 9 2 (1 - C 2 ). 



(57) 



(58) 



(59) 



Now, from the expression of L (56), we can also write 

,2 



L = w — (w 
Finally, we have 



and 



1/ = 



C 2 - 



C7(W 2 ) 
i-C 2 



(60) 



(61) 



and thus 



C= tvt(7 + fix + ay) 

Wq 



■ v / ap (3q 
Wo \ 7 7 



We have also from ( 33 1 



1 



'v 2 = ~2 (P 2 + q 2 + tyqz - (px + qy) 2 ) 



7 



so (61) can be written on the form 
9 2 = 0(x,y,z). 



(62) 
(63) 

(64) 

(65) 



The sign of 9 can be determined through ( 58 ) . Indeed 9 



is a function of (w, wj w? ), b ut its sign can only change 
when 9 = 0, that is from (61 1, when 



w 2 (i-C 2 ) = C 2 



(66) 



The equation ( 65 1 thus gives the instantaneous precession 



rate of w as a function of x, y, z. Same results can easily 
be obtained for the other two vectors Wi and W2. 

3.3.2 Symmetry of the nutation 

It is now possible to prove a more precise result on the 
periodic loops generated by w, W] and w 2 in the process- 
ing frame. This is the same result as in BL06 that was 
given for a three body problem with only one rigid body. 

Proposition 2 In the frame rotating uniformly with the 
precession period, the three vectors w, wi, W2 describe pe- 
riodic loops C, Ci, C2 that are all symmetric with respect 
to the same plane S containing Wo- 

Consequence. Let us call V, V\, V2 the averages of w, 
wi, W2 over the nutation angle. V, V\, Vi are respec- 
tively the pole of the orbit, the pole of the spin of the 
primary and the pole of the spin of the secondary. Due 
to the symmetry of the loops, the three poles V, V\ and 
V2 remain in the symmetry plane S containing wo, and 
precessing uniformly around wq. Each vector w, wx, w 2 
nutates around its pole, respectively V, Vx, Vi- 

Proof. As in BL06, we will consider uniquely w, the 
other cases being similar. We consider here a general so- 
lution, for which the orbit of (x, y, z) crosses the Cassini 
berlingot (Fig. |5|. We choose the origin of time in r + 
which corresponds to an orbital angular momentum w + . 
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Let a be the arc length described by M = (x, y, z) com- 
puted from M + = M(t + ). From (36 1 we have 



v\/ f(x,y,z) 



where 

f(x,y,z) 



(67) 



(68) 



f(x, y, z) = if and only if ap — [3q — js. This condition 
corresponds to a fixed point of the system. Else f(x 7 y, z) 
is strictly positive. 

Thus a is a function of (x,y 1 z) and has the sign of 
v. For t < 0, the orbit in the (x, y, z) describes the arc 
(t_,t + ), thus a decreases from er_ down to er + = 0, and 
v < 0. Conversely, for t > the orbit describes the same 
arc in the reverse way (t + ,t_), hence v > 0. As x, y, z 
are functions of the arc length a, we can write 



if * <0, 
if t > 



(69) 



where F(a) = \v\^J fix, y, z). We conclude that a and 
thus M = (x,y, z) are even, that is M(-t) = M(t). 
The rest of the proof is identical to the one of BL06. 



We recall it for completeness. From (651 
2 (t) = e(x,y,z), 



(70) 



we deduce that 9 2 (t) is even. Moreover, as the differential 
system (33 1 is polynomial, the solutions w, w 1; w 2 are 
analytical in time t, and so will be the coordinate angle 
9(t) of w. The lemma of BL06 thus implies that 9(t) is odd 
or even. If 9(t) is even on [-T/2, T/2], for all h € [0, T/2], 
we have 9(h)- 9(0) = 9(0)-9(-h). As the cosine C of the 
angle from w and w (55 1 depends only on x, y ( |62[ ), we 
have C(h) = £(— /i), and w(h) and w(— h) are symmetrical 
with respect to the (wo,w + ) plane. It will still be the 
same in the rotating frame with the precession period. In 
this rotating frame, the periodic loop generated by w is 
thus symmetric with respect to the plane (wo,w+). 

Moreover, at t = (r+), the volume v is null, and thus 
wo, w, wi, W2 are coplanar. In the rotating frame, all 
three orbits generated by w, Wi, w 2 are thus symmetrical 
with respect to the same plane (wo,w + ). 

The only case where 9(b) is odd, occurs when 9(0) = 0. 
As u(0) = 0, we have C(0) = {63) and w = 
In the same way, we have wi(0) = W2(0) = 0, and the 
vector field (33 1 vanishes at t = 0. The three vectors w, 
wi, W2 are thus stationary and coplanar. 

This is a special Cassini state where the precession fre- 
quency is zero. 



3.4 Computation of the two periods 

The nutation period and the precession period are two 
key parameters of the problem since the global solution 



is the product of these two motions ( 53 1 . Let us see how 
the values can be derived. 

The three dot products (x(t),y(t), z(t)) are T-periodic 
where T is the nutation period. This period can thus be 
calculated from the expression of (x(t),y(t), z(t)). Given 



the two first integrals (351, it is possible to express x(t), 
y(t), and z(t) in the form of an integral as in BL06. Nev- 
ertheless the energy conservation only gives an implicit 
relation between those variables and the computation re- 
mains tedious. For this reason, we give here an algorithm 
that enables to compute the two frequencies in a simple 
way using the numerical integration of the secular equa- 
tions ( 33 1 . The method leads to an arbitrary high preci- 



sion since it necessitates the integration over one nutation 
period only. 

We assume that at to = 0, the initial volume v (36 1 
is not zero, and let x (for example) be the variable with 
the largest variation rate, x(to). Using the method of 
Henon (1963), we search for the first time t > to when 
(x(t),x(t)) = (x(to),x(to)). We integrate the system (33 1 
until 



X n -l < x jf ± ( tQ j > o 

X n > X 



%n-l>Xa Hx(t )<0 
X n < Xq 



(71) 



(72) 



We then change the time variable to x and integrate 



dt 


1 


dx 


x(x,y,z)' 


dy _ 


y(x,y,z) 


dx 


x(x,y,z)' 


dz 


z(x,y,z) 


dx 


x(x,y,z)' 


d9 


(^J®(x,y,z) 


dx 


x(x,y,z) 



(73) 



from x n to Xq . The latter equation comes from ( 65 ) and 



will provide the rotation angle of the vectors over one 
nutation period (knowing the initial angle 9 (to)). We thus 
have the nutation period t = T and 9t = 9(T) — 6(to)- 
The precession period is simply given by 



r 



2tt, 



-T . 



(74) 
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4 Analytical approximation 

In this section we give an analytical approximation of the 
secular evolution. So far, only general features of the solu- 
tions have been obtained. Here analytical approximations 
of the two frequencies that appear in the problem as well 
as their amplitudes are computed. The two frequencies 
being the global precession and the nutation. 

In an invariant frame where the third axis is aligned 
with the direction Wo of the total angular momentum, we 
can write 




Wi 



w 2 



(75) 



where 



c = 



7 + (3x + ay 
^W> ' 



C-2 



Wo 



(76) 



The evolution of the projections on the complex plane 
orthogonal to wq 



3 = £ + 3X = £x + %, 32 = & + W72, 



(77) 



is obtained from the secular equations ( 33 1 , and yields to 

(78) 




= iM 



where 



M 



■fCx 



0^ 

*6 



7^ 

■|C-fC 2 



7 s 



;C- fCx 



(79) 



and (p, q, s) are defined in (34 1. M is a real matrix with 
periodic coefficients. As it is not possible to obtain a 
simple analytical solution of this system, we make a crude 
approximation. Hereafter we replace the matrix M by the 
constant matrix M obtained by substituting (x, y, z) by 
their average 

M = M(x,y,z) . (80) 

The solution of ( [78] ) is thus straightforward. It is easy to 
verify that (C, C11C2) is an eigenvector of M with eigen- 
value 0. The other eigenvalues are then the solutions of 



A 2 - TA + P 
where T is the trace of M and 

' C , Cx , C2 







p = 



s a(3 7a /?7 



(pq( + spd + qs( 2 ) 



(81) 



(82) 



Let and £1 + v be the other two eigenvalues such that 



T + VT 2 - 4P / 

0= 2 — , 1/= - a/T2 -4P . (83) 

The system possesses three eigenmodes 

ue^, te l ( nt+ *\ Se iW+u)t+$+< t >} ) ( 84 ) 

with eigenvectors 




where A, A', /t and // are real numbers. The solutions are 
then 

3 = Cue** + e^ t+ *> (r + sS vt+ ^) , 

31 = Cxue^ + e*( at+ *> (At + A'sc^ t+ «) , (g6) 

32 = C 2 uc 1 ^ + S nt+ ^ (/it + /i'se 4 ^*+«) . 



Moreover, 73 + (3$i + a$2 — as it is the projection of 
Wo on a plane orthogonal to V^o- This implies that 
the constant term (7^ + (3(i + aC2)ue 1 ^ is also null. As 
7C + (3(i + a( 2 — Wq, we have necessarily u = 0. The 
solutions are thus 



3 =c l ( ot +*)fr + se 1 ^ 



31 =e i{nf+*)J Ar+ y Se i(.i+flj ] 

32 = eW*)( M t + M '5e^ t+ ^). 



(87) 



In this approximation, the three axes (w, wi, W2) describe 
circular motions with nutation frequency v around the 
three poles (V,Vi,V2) that precess uniformly with pre- 
cession frequency f2 around the total angular momentum 
Wq. As it was previously said, the three poles (V, V\, V2) 
remain always coplanar with Wq. 

4.1 Initial conditions 



The preceding section shows that the solutions (87 1 de 



pend only on four real numbers r, s, $ and <f>. At the 
origin of time (t = 0) we can choose two vectors, for in- 
stance 



3o 



'*(t + and 310 = e l * (Xt + A'se 1 



from which we derive 



A' -A 



4, _ A30 - 310 
A -A' ' 



(89) 
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The computation of A and A' requires the knowledge of 
the averaged values x, y and z, but it can easily be done 
by iteration, starting with the initial values, that is, for 
the first iteration 

x = x(t = Q), y = y(t = Q), z = z(f = Q). (90) 

In our computations, we found that one iteration after this 
first try with the initial conditions was sufficient to obtain 
a satisfactory approximation for the frequency amplitudes 
and phases of the solution (see Tables [4] [7] [8]) . 

4.2 Second order expansion 

The whole previous study has been made with an Hamil- 
tonian expanded up to the fourth degree in R/r (TtT) , 
(pi p| and p} 



H = H T + H E + H 



(o) 



H 



(2) 



H 



(4) 



(91) 



When the body-body interactions are neglected, we can 
restrict the analysis to the second degree in R/r. The 
secular Hamiltonian then simplifies to 



— s 2 



Cte 



where 



and with 



a = k^m2C' 1 and b = hj,m\C' 2 
= 3 g 

2a 3(l- e 2)3/2 



(92) 

(93) 



C 1 =(l-| B iii 9 J 1 



Co 



1 



3 • 2 j 



C 2 



Ai + Bx 
2 

A 2 + B 2 



(94) 



The secular equations ( 33 1 become 

by 



ax 

w = — — Wl x w — 

7 



7 



-W2 X W 



Wl 



w 2 = 



ax 

'1 



W X Wl 



(95) 



-W X W2 



where 7, /3 and a are still the angular momentum of the 
orbit, of the rotation of the primary and of the rotation 
of the secondary respectively. In that case, the matrix M 
giving the evolution of the projection of the three vectors 
3, 31 and 32 becomes 



1/ = i f d -f C 

#C2 -& C y 



(96) 



Now we use the same trick as in the equation ( 80 1 , that 



is we replace the matrix M by the constant matrix M 

M = M(x,y,z) (97) 

where (x, y, z) have been substituted by their average. 
The vector *(£, (1, £2) is still an eigenvector for the eigen- 
value 0. The characteristic equation is now 



A - TA 







where 



ax by ax by 

S.1 <,2 C, 

7 7 p a 



P = abxyC 



C , Ci , C2 



(98) 



(99) 



7a /37 y 

These expressions give simpler formulas for the frequen- 
cies, although they still have the same form 



T+ VX 2 - 4P 



4P 



(100) 



5 Global precession of a n-body 
system 

We have seen that the secular motion of a two solid body 
system can, as in BL06, be decomposed in a uniform pre- 
cession of angular motion f2, and a periodic motion of fre- 
quency v. In fact, this can be extended to a very general 
system of n solid bodies in gravitational interaction. The 
following result, which is of very broad application, is a 
consequence of the general angular momentum reduction 
in case of regular, quasiperiodic, motion. 

Proposition 3 Let S be a system of n + 1 bodies of 
mass rrii,(i = 0, ...n) in gravitational interaction, with 
n s solid bodies among them (n s < n + 1). Then, in 
a reference frame centered on one of the bodies, and 
for a regular quasiperiodic solution of S , there exist a 
constant precession rate £1, such that any vector Z G 
{Ti,Ti,Ij,Jj,Kj,Gj]i = l,...n;j = l,...n s } has a 
temporal evolution that can be decomposed as 



z(t) = n 3 (nt)z^(t) 



(101) 



where 7^.3 (Ot) is a uniform precession around the total 
angular momentum Wq with constant rate fl, and where 
Z^\t) can be expressed in term of quasiperiodic series of 
3(n + n s ) — 2 frequencies (ffe). We will call £1 the global 
precession rate of the system S. 

Proof. Let us consider a general system of n + 1 bod- 
ies of mass rrii, (i = 0, . . . n) in gravitational interaction, 
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with n s solid bodies among them (n s < n + 1). This is a 
3(n + 1 + n s ) degree of freedom (DOF) system. Due to 
the translation invariance of the system, it can be reduced 
to N = 3(n + n s ) DOF using the coordinates centered on 
one of the bodies (the one of mass mo for example) . This 
heliocentric reduction can be made in canonical form, pre- 
serving the Hamiltonian structure of the equations (see 
Laskar and Robutel, 1995). 

The full Hamiltonian of the system, as ex- 
pressed in is then a function of the vectors 



l,...n;j 



1. 



that 



depends uniquely of the scalar products of theses vectors. 
Moreover, the total angular momentum Wo (35 1 is 
conserved. 

This system, as for the usual reduction of the node, can 
be reduced to a system of N— 2 degrees of freedom. A first 
reduction to N — 1 DOF can be achieved by using a refer- 
ence frame (i,j,k) such that k is collinear with Wq and 
k- Wq is positive. This partial reduction is based uniquely 
on the fixed direction of the angular momentum (Malige 
et al., 2002). With this reference frame, all quasiperiodic 
solutions of the system can be expressed in term of only 
N — 1 fundamental frequencies. 

In this fixed (i, j, k) reference frame, we can use canon- 
ical coordinates that are well adapted for both the orbital 
and rotational motions. Namely, we shall use the Andoyer 
coordinates for the solid bodies (L, G, H, l,g, h) (Fig. [2|, 
and the equivalent Delaunay coordinates for the orbital 
motions (A 



(3y/JI5,T = A\/l - e 2 ,0 = rcosi,M,a;,0) 
where (a,e,i,M,ui,9) are the usual elliptical elements 
(semi-major axis, eccentricity, inclination of the orbit with 
respect to the plane, mean anomaly, argument of pe- 
riapse, longitude of the ascending node). For any given 
body of mass m,,z ^ , 0i — momi/(mo + mi) is the 
reduced mass, and /ij = G(mo + mi) the related gravi- 
tational constant. For any Xi G {rj,fi;i = l,...n}, or 



Yj e {I j ,J j ,K j ,G j ;j=l, 



5 }, one can then write 



X i = n 3 {e i )X' i (A h T i> Q i ,M i ,uj i ) 
Yj =n 3 (h j )Y!(L j ,G j ,H j ,l j , 9j ) 



(102) 



Let us now select one angle among the 9i,hj (9\ for 
example) and perform the usual symplectic linear change 
of variable 



h'j = hj - 0i 



i 3 

©J: = © 4 for i ^ 1 
H'j = Hj 



(103) 



As the Hamiltonian (|T|) depends only on the scalar 
products of Xi and Yj, it can be as well expressed in 



term of scalar products of 



Yj = n 3 {-ei)Yj 



(104) 



Expressed in term of the new variables (1031, one can 



see that the coordinate 9[ is now ignorable with an as- 
sociated constant action being the modulus of the to- 
tal angular momentum (0^ = ||Wo||). The number of 
DOF of the system, expressed in the new coordinates 
(Ai^Q'^M^uj^e^Lj^j^'jJj^j^'j) is now N - 2, 
with one constant parameter, 0' x . Let us now consider 
a quasiperiodic solution of the above N — 2 DOF sys- 
tem. All vectors Xi,Yj will be expressed in term of 
quasiperiodic functions on N — 2 independent frequen- 
cies i>k, {k = 1, ■ ■ ■ N — 2). Finally, 9[ evolution is given 

by 



d6[ _ dH 



(105) 

Thus 9[ (t) is also a quasiperiodic expression depending 
on the N — 2 frequencies Vk . 



d9[ 
~dt 



exp(i < k,v > t) , 



(106) 



where (fc) is a (N — 2) multi index. Let f2 = a( ) be the 
constant term of this series. We have then 



d9[ 
and thus 



= Q+ V" exp(i < fc, v > t) , (107) 

(fe)^(O) 

6[(t) =m + f (v) (t) , (108) 

where f( v )(f) is a (N — 2)— periodic function with fre- 
quencies (ffc). The original vectors X i: Yj can then be 
expressed as 



W 

i 

O) 



Xi - n^Xi = n 5 {nt)n 3 {f {u) {t))Xi = n 3 (nt)x 

Yj =7£ 3 (0i)^ =K 3 (nt)1l 3 (f {u) {t))Yj = K 3 (nt)Y 

(109) 

where x\ v \y^ v) can be expressed in term of (N — 
2)— periodic function with frequencies (vk)- This ends the 
proof of the proposition. 

Consequence. A consequence of this result is that for 
a quasiperiodic solution of the general two body problem 
that we are considering here (n — l,n s — 2), the com- 
ponents of any vectors r,r,Ij,Jj, Kj , Gj , should express 
as quasiperiodic functions of the precessing frequency fl 
and of 7 frequencies Vf.,k = 1, ...7, the precession fre- 
quency f2 appearing in all terms with coefficient 1. This 
is actually what is observed on some examples in the next 
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section (Tables [5] and [6]) . One should note that the same 
results hold for the three body problem studied in BL06 
(with n = 2, n s = 1). 

It is also useful to remark that the value of is indepen- 
dent of the i>]~, i.e. any commensurable relation between 
f2 and the Vk has no effect on the dynamics of the system, 
in the sense that it will not affect the regularity of the 
solutions. On the other hand, in the case of a single Vk 
frequency (as for the secular system) , a rational ratio Cljv 
will lead to a periodic solution in the fixed reference frame 
(i,j,k). We prefer here to speak of geometric resonance 
instead of dynamical resonance, as there is no coupling 
between the two degrees of freedom of frequency f2 and 
v. 



6 Application 

In this section we compare our rigorous results on the av- 
eraged system and our analytical approximations of the 
solutions of the same system with the integration of the 
full Hamiltonian (|2]), 0, ([l3]) and on two dif- 

ferent binary systems I and II (see table^TJand [2]) . The 
physical and orbital parameters of the system 17 are those 
of the binary asteroid 1999 KW4 studied in FS08. We 
choose this system in order to compare our results with 
FS08. In this case, the rotation of the satellite is taken to 
be synchronous. As our analytical results were obtained 
assuming the satellite rotation asynchronous, we create a 
system / from the system II where the rotation of the 
secondary has been sped up by a factor 3. Since the or- 
bit is circular and the initial rotation axes aligned with 
the axes of maximum inertia, the system II is highly de- 
generated. To get a more general system where all the 
fundamental frequencies will actually exist, we changed 
the initial Andoyer angles and the eccentricity. But then, 
because of its strong triaxiality, the evolution of the satel- 
lite orientation becomes chaotic (Wisdom, 1987). As here, 
we are concerned only with on regular behaviors, we thus 
decreased the satellite triaxiality and increased the semi- 
major axis in order to obtain a generic example of regular 
solution. 



6.1 Numerical experiments 
6.1.1 Frequency analysis 

The quasiperiodic decomposition of our numerical inte- 
grations was obtained using the frequency analysis devel- 
oped by Laskar (Laskar, 1988, 2005). As our systems con- 
tain a large range of frequencies going from 0.07 rad-day -1 
to 109 rad-day -1 , we decided to run twice each integra- 
tion with two different output time steps ft, = 0.1 days and 



Table 1: Physical and orbital parameters of a fictitious doubly asyn- 
chronous binary system, m is the mass, A, B and C are the moments 
of inertia divided the mass, w is the rotation rate, h, 7, g, J and I 
are the Andoyer angles of the two solid bodies as defined in Fig. [2] 

System I 





Primary 


Secondary 


Orbit 


m (10 12 kg) 


2.5 


0.15 


a (km) 


2.75 


A (km ) 


0.17 


0.0165 


A (deg) 


0.0 


B (km 2 ) 


0.18 


0.017 


e 


0.035 


C (km 2 ) 


0.19 


0.025 


w (deg) 


0.0 


w (°/day) 


3125.34 


1500 


i (deg) 


0.0 


h (deg) 


100.82 


-110.0 


d (deg) 


0.0 


/ (deg) 


10.74 


5.0 






9 (deg) 


112.03 


-180.0 






J (deg) 


3.0 


5.0 






I (deg) 


90.0 


90.0 







h' = 0.1001 days. These two time steps do not fulfilled 
the Nyquist condition for the largest frequency. Never- 
theless, it is possible to recover the true value v$ of the 
frequency using the following trick (Laskar, 2005). For a 
real x, let denote [x] the real such that 



7T < [x] < 7T. 



(110) 



Let v and v' be respectively the frequencies measured on 
the integration with the time step h and h! . The true 
frequency is given by 



h 



where 



-{(v' -v)ti - [v'h'-vh]). 



(Ill) 



(112) 



6.1.2 System I doubly asynchronous case 

Full Hamiltonian We integrated the system I over a 
time span of 2 000 days and performed a frequency anal- 
ysis as described above. This system contains a priori 9 
degrees of freedom. Three coordinates for the orientation 
of each body and three coordinates for the orbit. But 
there is a relation between all these coordinates given by 
the conservation of the total angular momentum. There 
are thus only 8 degrees of freedom. Hence the system 
contains 8 fundamental frequencies (cf table [3]) . 

These frequencies can be divided into four main cate- 
gories: 1) the secular frequencies containing the preces- 
sion CI and the nutation v\ 2) the orbital frequencies with 
the periapse precession rate uj and the mean motion n; 
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Table 5: Frequency decomposition of the motion of the projections 3, 31 and 32 respectively of w, wi and W2 on the plane orthogonal to the 
total angular momentum Wo- The integration was made using the full Hamiltonian with the initial conditions of the doubly asynchronous 
system I. Only the first 20 terms of the series Aj expi(i>jt + ipj) are displayed for each vector. In order to simplify the reading, hats on 
the angles cj, gi, li, 32 and I2 are omitted. 



var. 






/_ „ j ,„,-ii 
(lad.yi _) 


( ) 


(dcg) 




1 


1 1 


-0.0312 


on 1 n k nn 
Z\) lOo.Oy 


-169.37 




2 


Q + v 


-1.0100 


209.11 


-17.70 




3 


fl + 2uj + 2n 


16.1156 


55.87 


-10.55 




4 


Q — v + 2lj + 2n 


17.0943 


12.29 


-162.22 







ii zw -j- zn — p2 


-zo.oo4o 


9.55 


-158.86 




6 


ii -j- 2cj + 2n — gi 


-42. 8607 


6.70 


57.27 




7 


f2 + p2 


39.9391 


5.72 


-21.06 




8 


Q — 71 


-8.0365 


5.22 


-168.69 




Q 


Q _|_ n 


7 9740 


5 21 


9 94 


w 


10 


n + pi 


CO n A K 1 
o8.94ol 


4.89 


122.82 


1 1 


n + 2cj + 3n 


24.1208 


4.00 


-11.23 




12 


ii -j- Jo; -+- zn — zii — zpi 


-y3.024o 


3.39 


-55.00 




13 


SI + 2/i + 2pi 


1 nn inon 


2.90 


-124.92 




14 


ii -+- -+- n 


8.1103 


1.70 


170. 13 




15 


li zw -J- ofi — p2 


-15.8495 


1.47 


-159.55 




16 


li -+- zlj -j- Zn — £l\ — Q\ 


-34.0400 


1.39 


57.17 




17 


Sa + v 4- n 


6.9953 


1.08 


161.62 




18 


Vl -\- V — 71 


-9.0152 


1 .04 


-17.02 




19 


s~2 — v 


0.9476 


0.95 


-141.05 




20 


S"2 4- 2l x 4- p x 


O0.132O 


0.94 


122.90 




1 




-0.0312 


9687.25 


10.63 




2 


n + 2oj + 2n 


16.1156 


18.80 


169.45 




3 


O + 2a; 4- 2n — pi 


-42.8607 


2.26 


-122.72 




4 


Q + v 


-1.0100 


1.86 


162.30 




5 


S7 — ri 


-8.0365 


1.73 


11.31 




6 


O 4- n 


7.9740 


1.73 


-170.06 




7 


ft + gi 


58.9451 


1.60 


-57.20 




8 


O + 2oj + 3n 


24.1208 


1.34 


168.77 






li — ZUJ — Zn -\- Zl\ -\- zg\ 


y z . ijoz 1 


111 


1 O.zO 


Wi 


10 


SI + 2l\ + 2pi 


1 nn inon 
lUy.lUoy 


0.94 


55.08 


1 1 


ii -f- zcj -f- n 


8.1103 


0.58 


-9.87 




12 


ii -+- Zuj -+- zn — Zl\ — g\ 


OA n A Q "3 

-34.U483 


0.47 


-122.83 




13 


Q + 2/i + pi 


O0.13zo 


0.31 


-57. 10 




14 


li -f- ZCJ -f- OTI — pi 


-34.8000 


0.30 


- 123.41 




15 


O — 2a; — 2n 


1 a 1 Ton 
-16.1 r oO 


0.16 


- 148. 19 




10 


ii — zo; — on -\- zi\ -+- zg\ 


o4.yooy 


0. 13 


76.94 




17 


ii — y + ^oj + zn 


1 T nn /i "3 
1 I .0y43 


0.13 


17.78 




18 


f2 — n + pi 


00.939o 


0.08 


-56.52 




10 


c"^ 1 r> H , 1 j „ 
1 i + zlj + 4n 


go 1 ocn 


0.07 


168.08 




20 


ii + zlj -j- on — zti — pi 


oc n /i q n 
-20.0430 


0.07 


- 123.51 




1 




-0.0312 


■3 n n 7n 10 

30079.18 


-169.37 




2 


O 4- v 


-1.0100 


1 7701 00 
1 ( /81.88 


162.30 




3 


Q, — v + 2lj + 2n 


17.0943 


1040.78 


17.78 




4 


SI + 2uj 4- 2n — p2 


-23.8548 


829.13 


21. 13 




5 


H + p 2 


39.9391 


495.19 


158.94 







O + 2cj 4- 3n — p2 


-15.8495 


126.95 


20.45 




7 


Q — 1/ 


0.9476 


96.54 


38.96 




8 


i2 -\- v -\- 71 


6.9953 


90.98 


-18.39 




Q 


1 i -r f — J t 


-9 0152 


87 74 


162 98 


w 2 


10 


S7 — v 4- 2w + 3n 


25.0996 


74.38 


17.10 


11 


S"2 — f 4- 2w + 71 


9.0891 


50.41 


-161.54 




12 


n + 2cj + 2?i 


16.1156 


46.84 


-10.55 




13 


£1 — n + p2 


31.9339 


25.33 


159.63 




14 


tt + 2w + 2n - 2l 2 - P2 


3.9537 


23.66 


-161.23 




15 


Q + 7i 4- p 2 


47.9444 


21.72 


158.26 




10 


S"2 4- 2a; + 4n — p 2 


-7.8443 


19.43 


19.77 




17 


SI + 2n + 2/i - 2/ 2 


34.9754 


14.90 


-57.53 




18 


n - y + 2/ 2 + 2p 2 


53.0797 


10.87 


-22.05 




10 


S^ + 2/ 2 + 92 


12.1306 


8.74 


161.31 




20 


Q, + 2w 4- 5n - p 2 


0.1610 


7.16 


-160.92 



3) the frequencies of the primary g\ and l\ associated re- 
spectively to the Andoyer angles g\ and 1%; 4) the same 
frequencies for the secondary §2 and li- 



Table [5] displays the frequency decomposition in the 
form ^2 Aj exp i(vjt + ipj) of the motion of 3, 31 and 32 
(77), the projections of w, wi and W2 on the complex 
plane orthogonal to W . The second column shows that 
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Table 6: Same as table 5 for the system II. 



var. 


i 






v i 


(") 


<p\ f) 








(rad.yr- 1 ) 


(dog) 




1 


fi 


-0 


0713 


27979.553 


-153.15 




2 


SI + 2u + 2n 


17 


8490 


111.214 


-26.81 




3 




_4 


8201 


4.421 


-160.28 




4 


Q + 2ui + 2n - 2h - 2gi 


-91 


3906 


3.320 


-39.33 




5 


SI — n 


_g 


1216 


2.847 


25.69 




6 


S2 + 2ii + 2gi 


109 


1682 


2.779 


-140.63 




7 


n + n 


g 


9790 


2.678 


-151.99 




8 


SI + 2ui + 3n 


20 


8993 


2.201 


154.37 




9 


fi+ V>2 




5201 


1.295 


-152.20 


w 


10 


n- v>2 


-7 


6627 


1.117 


25.96 


11 


Q + 2w + n 




7986 


1.066 


-27.97 




12 


SI - i/ + 2w + 2ri 


22 


5978 


0.942 


-19.67 




13 


SI + 2a; + 2n + i/>2 


25 


4404 


0.828 


155.16 




14 


n - 2a! - 2n 


-17 


9916 


0.490 


-99.49 




IS 


S2 + oj + n — 0-2 


4 


7413 


0.376 


-163.33 




16 


n + 2lo + 2n - lp 2 


10 


2575 


0.280 


-27.75 




17 


n + lo - 2 


_4 


3090 


0.237 


-164.49 




18 


S2 — v + 2lj — n 


_4 


5532 


0.170 


156.75 




19 


SI - v + 2lo 


4 


4971 


0.161 


157.74 




20 


S2 + f + n 


4 


2302 


0.143 


20.77 




1 


SI 


-0 


0713 


8008.982 


26.85 




2 


Si + 2u + 2n 


17 


8490 


32.172 


153.19 




3 


n - 2uj - 2n + 2l x + 2g t 


91 


2480 


0.939 


93.03 




4 


VL — n 


_g 


1216 


0.796 


-154.31 




5 


SI + n 


g 


9790 


0.794 


28.01 




6 


S2 + 2fi + 2gi 


109 


1682 


0.780 


39.37 




7 


U + 2u + 3ti 


26 


8993 


0.636 


-25.64 




8 


n - ip 2 


_7 


6627 


0.325 


-154.10 




9 


SI + l/>2 


7 


5201 


0.324 


27.80 


1 


10 


SI + 2w + n 


8 


7986 


0.301 


152.03 


11 


S2 - 2w - 2n 


-17 


9916 


0.243 


-99.49 




12 


SI + 2a; + 2n + i/>2 


25 


4404 


0.232 


-25.86 




13 


S2 + 2w + 2n - V> 2 


10 


2575 


0.093 


152.25 




14 


S2 - 2w - 3n + 2ii + 2gi 


82 


1976 


0.031 


-88.11 




15 


SI + ^ 


_4 


8201 


0.030 


19.72 




16 


S2 + 4w + 4n - 2ii - 2gi 


-73 


4703 


0.029 


87.00 




17 


SI + 2w 


-0 


2517 


0.024 


-29.13 




18 


S2 + 2lo + 2ri - 22 x - 2g t 


-91 


3906 


0.023 


140.67 




19 


il + n — tp2 


1 


3876 


0.020 


-152.96 




20 


£1 — n + ip 2 


-1 


5302 


0.020 


26.63 




1 


n 


-0 


0713 


28848.685 


-153.15 




2 


S2 + v 


_4 


8201 


924.226 


19.72 




3 


Q - v + 2io + 2n 


22 


5978 


196.275 


160.33 




4 


SI + 2u) + 2n 


17 


8490 


81.967 


-26.85 




5 


SI + u + n — 82 


4 


7413 


77.510 


16.67 




6 


S2 + v + n 


4 


2302 


55.301 


-159.15 




7 


S2 + a; - e 2 


_4 


3090 


47.469 


15.51 




8 


S2 — ^ + 2^ — n 


_4 


5532 


34.386 


-23.16 




9 


U - v + 2lu 


4 


4971 


33.518 


-22.15 




10 


SI + ^ — 1p 2 


-12 


4116 


30.128 


18.77 


Wo 


11 


n + ui + n + 9 2 


13 


0363 


29.773 


-16.63 




12 


n - -!/> 2 


-7 


6627 


29.524 


-151.52 




13 


S2 + w + 6 2 


3 


9860 


26.699 


162.21 




14 


Q -\- v — n 


-13 


8705 


26.395 


-161.45 




15 


il + u - n + 2 


-5 


0643 


19.869 


-18.95 




16 


a + v + v> 2 


2 


7713 


19.196 


20.64 




17 


S2 + a; + n — i/>2 + #2 


5 


4449 


18.432 


162.43 




18 


S2 + 1/ + 2n 


13 


2805 


13.450 


22.33 




19 


U - v + 2iu + 2n - V>2 


15 


0063 


12.645 


159.39 




20 


Q. + LO + U — \p2 — 02 


-2 


8501 


12.132 


-164.27 



all the frequencies are combinations of the 8 fundamental 
frequencies. 

Moreover, we verify our proposition saying that in a 
frame rotating uniformly with the precession rate f2, the 
system loses one degree of freedom, see section [5] Indeed, 
the frequency £1 appears in all the terms with the same 



order 1. 



Averaged Hamiltonian In the frequency decomposi- 
tion of the motion of Wi in Table [5j the nutation is only 
the 4 th term. To check the validity of the averagings, 



we integrated the averaged Hamiltonian ( 28 1 on the same 
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Table 7: Frequency analysis of the doubly asynchr onou s system /. Columns 3 to 5 correspond to the frequency analysis performed on the 
numerical integration of the averaged hamiltonian |28| , superscript (a). Columns 6 to 8 contain the secular terms of the frequency decom- 
positions computed on the output of the full integration, superscript (/). Columns 9 to 11 are the results of the analytical approximations 
l |83| and |89| , superscript (c). 



var. 






i 


( ) 


a 

(deg) 


% 


(II \ 

( ) 


(/) 
¥>< 

(deg) 


i 


Ac) 

fii\ 

( ) 


(c 

(deg) 






9 


1 


29104.12 


-169.37 


1 


29105.09 


-169.37 


1 


29108.16 


-169.37 


w 


9 


+ v 


2 


209.34 


-17.37 


2 


209.11 


-17.70 


2 


212.43 


-17.10 


fl 


— V 


3 


0.95 


-141.38 


19 


0.95 


-141.05 










9 4 


- 2v 


4 


0.04 


134.64 


57 


0.04 


134.09 












9 


1 


9688.12 


10.63 


1 


9687.25 


10.63 


1 


9688.16 


10.63 


Wl 


n 


+ v 


2 


1.76 


162.63 


4 


1.86 


162.30 


2 


1.85 


162.90 




9 


— V 


3 


0.05 


-141.38 


23 


0.05 


-141.02 












9 


1 


30020.39 


-169.37 


1 


30079.18 


-169.37 


1 


29676.07 


-169.37 




9 


+ v 


2 


17889.69 


162.63 


2 


17781.88 


162.30 


2 


18137.11 


162.90 


w 2 


9 


— V 


3 


96.21 


38.63 


7 


96.54 


38.96 










fM 


-2v 


4 


3.29 


-45.36 


37 


2.95 


-46.02 










9 - 


-2v 


5 


0.02 


-113.38 















Table 8: Same as tabic 7 for the system II corresponding to the 1999 KW4 binary asteroids. 



var. 






i 


(") 


(a) 

(deg) 


i 


A U) 
(") 


(deg) 


i 


(") 


(deg) 






9 


1 


27916.13 


-152.96 


1 


27979.55 


-153.15 


1 


27916.04 


-152.96 


w 


n 

9 


+ v 

— V 


2 
3 


4.65 
0.02 


-152.96 
27.04 


3 


4.42 


-160.28 


2 


4.72 


-152.96 


Wl 




9 


1 


7991.22 


27.04 


1 


8008.98 


26.85 


1 


7991.22 


27.04 


9 


+ v 


2 


0.04 


27.04 


15 


0.03 


19.72 


2 


0.04 


27.04 






9 


1 


28903.02 


-152.96 


1 


28848.69 


-153.15 


1 


28922.58 


-152.96 


w 2 


9 


+ v 


2 


987.12 


27.04 


2 


924.23 


18.39 


2 


1001.81 


27.04 


9 


— V 

-2v 


3 

4 


4.87 
0.01 


-152.96 
27.04 


53 


1.48 


-146.03 









time span (2 000 days) and we performed the same fre- 
quency analysis. Initial rotation rates, semi-major axis 
and eccentricity are average values computed on the nu- 
merical output of the full integration. Initial inclination, 
obliquities and ascending nodes were obtained from the 
amplitudes and the phases of the frequency analysis in 
Table E 

Table [7] displays the comparison between the frequency 
decomposition of the output of the averaged Hamiltonian 
and of the full Hamiltonian (columns 3-8). For the com- 
parison, only the secular terms were extracted from the 
analysis of the full integration. The second column con- 
firms our analytical result saying that the averaged mo- 
tion contains only 2 fundamental frequencies: the preces- 
sion 9 and the nutation v; and that in a frame rotating 
with the precession frequency, only the nutation remains. 



The columns 4,5 and 7,8 show the strong agreement be- 
tween the secular approach and the full integration. Even 
low amplitude terms such as 9 + 2v, albeit at the 57 th 
position in the decomposition of w in the full integration, 
are recovered with good amplitude and phase in the reg- 
ular system. 

The two last columns of table [7]give the complex ampli- 
tudes of the secular motion obtained with the analytical 
approximation of section [4] As in this approximation, the 
nutation is assumed to be a uniform rotation, there are 
only two terms in the description of the secular motion. 
Nevertheless, we see that this approximation is also in 
good agreement with the integration of the full Hamilto- 
nian and of the averaged Hamiltonian. 

In Tableware given the values of the secular frequencies 
for systems / and If, obtained either from the integration 
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x quasi-component (deg.) x quasi-component (deg.) 



Figure 6: Quasi-projection of the poles w (red), wi (green), W2 (blue) on the plane perpendicular to the total angular momentum V^Oi 
in a fixed reference frame (left panel) and in a frame rotating with the precession period (right panels). The two little figures on the right 
are zooms on the nutation motion of the orbit (top) and of the primary axis (bottom). The initial conditions and parameters are those of 
the system /. The vectors w, wi and W2 have been integrated with the full Hamiltonian. In the right panels, the output of the averaged 
Hamiltonian has been superposed: w in cyan, wi in pink and W2 in orange. 



Table 2: Physical and orbital parameters of the binary asteroids 
1999 KW4 given by FS08. m is the mass, A, B and C are the 
moments of inertia divided the mass, w is the rotation rate, h, I, g, 
J and I are the Andoyer angles of the two solid bodies as defined in 
Fig. g 

System II 





Primary 


Secondary 


Orbit 


777 (10 12 kg) 


2.353 


0.135 


a (km) 


2.5405 


A (km 2 ) 


0.1648 


0.01608 


A (deg) 


0.0 


B (km 2 ) 


0.1726 


0.02374 


e 


0.0 


C (km 2 ) 


0.1959 


0.02799 


ui (deg) 


0.0 


w (°/day) 


3125.34 


498.09 


i (deg) 


0.0 


h (deg) 


117.04 


0.0 


(deg) 


0.0 


/ (deg) 


10.0 


0.0 






9 (deg) 


0.0 


0.0 






J (deg) 


0.0 


0.0 






I (deg) 


-173.93 


180.0 







of the full Hamiltonian, from the integration of the aver- 
aged Hamiltonian, or with the analytical approximations 
(83 1. The precession rates arc in agreement within 0.3% 



and the nutation frequencies within 5%. 

Figure [6] represents the trajectories of the unit vectors 



Table 3: Fundamental frequencies of the two systems. Q and v 
are the precession and nutation frequencies respectively, ui and n 
correspond to the precession of the periastre and the mean motion. 
gi and ii on the one hand, and §2 and I2 on the other hand, are 
associated to the Andoyer angles. t/>2 and 82 are the horizontal and 
vertical libration frequencies in the resonant system i7. 

frequency value (rad/day) 

system / system II 



n 


-0.0312 


-0.0713 


V 


-0.9788 


-4.7488 


Co 


0.0681 


-0.0902 


n 


8.0052 


9.0503 


91 


58.9763 


63.3416 


h 


-4.4062 


-8.7218 


h 


39.9703 




h 


-13.9042 








7.5914 


02 




4.1475 
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Table 4: Secular frequencies. Comparison between the integration 
of the full Hamiltonian, the integration of the averaged Hamiltonian 
and the analytical approximations. 







n 


V 


system 


type 


(rad/day) 


(rad/day) 




full 


-0.0312 


-0.9788 


syst. I 


averaged 


-0.0310 


-1.1276 




calculation 


-0.0310 


-1.1091 




full 


-0.0713 


-4.7488 


syst. II 


averaged 


-0.0710 


-3.5982 




calculation 


-0.0710 


-3.5629 



on the plane orthogonal to the total angular momentum 
Wq. In the left panel, the frame is fixed, it thus corre- 



sponds to W(i), see section 3.2 We see that the evolu- 



tions of w in red and wi in green are dominated by the 
precession: their orbits are quasi-circular, whereas the or- 
bit of W2 in blue contains a large nutation as it can be 
checked in the frequency analysis table [5] The right panel 
shows the same orbits but in a frame rotating with the 
precession rate f2, it corresponds to W(t). It emphasizes 
the nutation loops. Zooms on the nutation of the orbit 
and of the primary are plotted on the furthest right. The 
solid curves are the output of the averaged Hamiltonian. 
The analytical approximations cannot be distinguished 
from the averaged output. These averaged solutions are 
good approximations of the motion of w and w 2 but the 
agreement does not seem to be as good for Wi. Indeed, 
table[5]shows that high frequencies have larger amplitudes 
than the secular nutation. 

Because of this amplitude issue for Wi in Fig. [6] we de- 
cided to filter our full integration with a low-pass filter to 
see if we could get back the averaged integration. In this 
scope, we reintegrated the full Hamiltonian over a time 
span of 20 days with an output time step of 30 min. We 
then filtered the output with a cutoff frequency equal to 4 
rad/day. The filtered trajectories are displayed in Fig. [7] 
The nutation amplitude of wi is now well retrieved. After 
a small change in the initial conditions that corresponds 
to a decrease of only 3.6" of the initial obliquity of the 
primary in the averaged Hamiltonian, we get back the 
filtered full Hamiltonian (see Fig. [8]). 



6.1.3 



System 
case 



II asynchronous-synchronous 



For this second experiment, we took the same initial con- 
ditions as FS08 (table [fjj. The primary has an asyn- 
chronous rotation whereas the secondary rotates syn- 
chronusly. The difference between our study and FS08 is 
that we expanded the Hamiltonian up to the fourth order 



in R/r where R is the radius of one body and r the dis- 
tance between them. We performed the same frequency 
analysis as with system /. We get also 8 fundamental 
frequencies. Because the resonance, the frequencies asso- 
ciated to the secondary are not <?2 and I2 anymore since 
they are in that case combinations of the other 6 funda- 
mental frequencies. The two new frequencies correspond 
to the horizontal and vertical libration of the secondary: 
1P2 and 82 respectively. 

Table [6] presents the frequency analysis performed on 
this system. The result of section [5] is still valid, in a 
frame rotating with the precession rate, the system loses 
one degree of freedom. We confirm that this result does 
not depend on the resonances in the reduced problem. 

The averaged Hamiltonian and the analytical approx- 
imation were not specifically written for such a resonant 
case. Regardless of this fact, the results of the averaged 
Hamiltonian and of the analytical approximation applied 
to this system are summarized in table [8] It is remark- 
able that the first two amplitudes of each vector are in 
good agreement with the full integration. Nevertheless, 
the third amplitude of W2 is wrong by a factor 3. The 
values of these secular frequencies are given at the bot- 
tom of table [4] The use of the averaged Hamiltonian or 
of the analytical approximations leads to an error on the 
precession rate equal to 1% and on the nutation rate equal 
to 24%. 

6.2 Comparison with FS08 
6.2.1 Numerical results 

In FS08, Fahnestock and Scheeres expanded the Hamil- 
tonian up to the second order in R/r. They find that 
motions of binary asteroids such as 1999 KW4 are combi- 
nations of four modes with their respective fundamental 
frequency. The first and fastest mode corresponds to the 
rotation of the primary around its axis. The second mode 
coincides with the orbital motion which has the same pe- 
riod as the rotation of the secondary around its axis. The 
third mode is said to be an excitation of the satellite's free 
precession dynamics and has a period of w 188 h. The 
corresponding frequency would be ~ 0.802 rad/day. The 
last mode is identified as the precession motion. 

Our results generally agree with the analysis of Fahne- 
stock and Scheeres. Nevertheless, several frequencies are 
missing in their analysis, probably because of the degen- 
eracy of their initial conditions. As the initial eccentricity 
is close to 0, and the angular momenta along the axes of 
maximum inertia, the first terms in the frequency decom- 
positions are combinations of Co + n which corresponds to 
their orbital frequency, and of g\ + l\ which corresponds 
to their rotation of the primary, see table [6] On the 
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0.08 
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-0.08 
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-2.6915 



-2.692 



-0.001 
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Figure 7: Same as the right panel of Fig. [S] The output of the full 
Hamiltonian, integrated over 20 days with an output step of 30 min, 
has been filtered with a low-pass filter with a cutoff frequency equal 
to 4 rad/day. 




-2.6915 - 



-2.692 



-0.001 







Figure 8: Same as Fig. [7] The initial obliquity of the 
the averaged Hamiltonian has been decreased by 3.6". 



with x = (w- Wi) and 



0.001 



primary in 



other hand, we do not find their third mode of frequency 
« 0.802 rad/day. 



6.2.2 Solid-point interaction 

Fahnestock and Scheeres found also that the spin axis of 
the primary and the orbital plane precess at the same rate. 
They derived an analytical expression for this precession 
rate, see their equation (76). Their result corresponds 
in fact to the solution of the single planet case that is 
already described in BL06 and which does not require the 
more elaborated formalism developed here. Indeed, as 
they expanded the potential up to the second order only, 
they canceled the effect of the orientation of the secondary 
on the precession of the primary (c = 5 = e = f = g = 0). 
Moreover, as they fixed the orientation of the secondary 
with the orbit, the secondary does not influence the orbit 
(y = (w- wj) = 1). We recall here the derivation of this 
frequency as given in BL06. With the assumption of a 
point mass satellite, the Hamiltonian becomes 



n 



(113) 



w = (w- w 1 )w 1 X w, 

7 



(114) 



Wl 



(w- Wi)w x wi. 



This reduced problem has 5 independent integrals given 

by 

1 

(115) 



w| 

|Wi| 



1 



7W + /3wi = Wq . 

As x = w • wi is constant, the system is trivially inte- 
grable. We have indeed 



w = fio w o x w i w i = ^o w o x wi 



(116) 



where Wo = Wo/ ||Wo|| is the unit vector in the direction 
of the total angular momentum Wq, and 



O 




(117) 



Both vectors w, Wi thus precess uniformly around the 
total angular momentum direction Wo with constant pre- 
cession rate f2o- The correspondence with the notations 
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of the equation (76) of FS08 is 



7 2a 7 /2(l_ e 2)2' 
x = cos((5 + i) 



(h -/)!-- sin 2 Ji 



(118) 




sin(<5 + i) 



sm z 



Remark. The factor (1 - (3/2) sin 2 Ji) is not in FS08 
because FS08 assumes that the primary angular momen- 
tum is aligned with its figure axis (Ji = 0). This is not 
the case in this paper where we do not require this sim- 
plification. 

7 Conclusions 

We have shown here that the general framework devel- 
oped in BL06 applies as well to the problem of two rigid 
bodies orbiting each other. This formalism enables us to 
obtain the long term evolution of the spin axis of the two 
bodies as well as the evolution of the orientation of the 
orbital plane. The two bodies can be very general, with 
strong triaxiality, and their rotation vector is not neces- 
sary aligned with their axis of maximum inertia. The 
gravitational potential is expanded up to the fourth order 
so as to keep the direct interaction between the orienta- 
tion of the two bodies, and as in BL06, the evolution of 
their spin axis is obtained after a suitable averaging. 

We found that the secular evolution is composed of two 
periodic motions: a global precession of the three angular 
momenta and nutation loops. As in BL06, the nutation 
loops are symmetric with respect to a plane containing the 
total angular momentum and precessing with the global 
precession frequency. We gave analytical approximations 
of these frequencies. 

We performed a frequency analysis (Laskar, 1988, 2005) 
on a numerical integration of the full Hamiltonian. We 
chose the typical binary asteroid system 1999 KW4 al- 
ready analyzed in FS08. We retrieved the precession and 
the nutation motions predicted by the secular Hamilto- 
nian and estimated by the analytical approximations. On 
a non resonant system, derived from 1999 KW4, the sec- 
ular solution, and the analytical results agree extremely 
well with the full solution. This is still the case to a 
lesser extent with the more specific case of 1999 KW4, 
which is in 1:1 spin-orbit resonance. In a further work, 
we could consider in a more precise way the possible res- 
onances. In that case, some of the averagings need to be 
done in a different way, probably leading to less symmet- 
rical, more complex, expressions. The main goal reached 



by the present paper was to search, in this apparently dif- 
ficult problem of two solid bodies in interaction, what was 
the most simple relevant underlying structure. One can 
now add possible additional effects, as tidal dissipation, 
and still consider the problem with the present setting. 
We thus expect that the results presented here will be 
helpful for the understanding of the general evolution of 
binary asteroids, or other problems of astronomical inter- 
est. 

In the elaboration of this paper, we came across the 
very general result given in our proposition [3] which ap- 
plies to any system of n massive bodies (point masses or 
not) in gravitational interaction. This property of the mo- 
tion states that the general regular quasiperiodic motions 
with N independent frequencies can be decomposed into 
a uniform rotation around the total angular momentum, 
which we call the global precession, and in this rotation 
frame, a quasiperiodic motion with N — 1 frequencies, in- 
dependent of the global precession frequency. 
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Appendix A. Gravitational interaction ex- 
pansion 



The gravitational interaction between the two bodies is 
given by (|5| 



Q dmi dm,2 
!|r + r 2 -nj 



(119) 
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The expansion of this potential in Legendre polynomials 
leads to the following integrals 



where X — x/a, Y = y/b, Z = z/c and B is the unit ball. 
From the definition of the moments of inertia 



H 
H 

H 



(i) _ 



1 

(2) 
I 

(3) 



= 



Q 

r 

Q_ 

2r 3 

Q_ 
2r 4 



Q 

8r 5 



dmi drri2 , 

3(u- ri) 2 + 3(u- r 2 ) 2 - r\ - r\ dmidm 2 
5(u-r 1 ) 3 -5(u-r 2 ) 3 
— 3r 2 (u- ri) + 3r 2 (u- r 2 ) dm,\dm,2 

35(u-r!) 4 + 35(u-r 2 ) 4 

-120(r r r 2 )(uT 1 )(u-r 2 ) 
-12r 2 (r 1 -r 2 ) + 210(u-r 2 ) 2 (u-r 1 ) 2 
-30r^(u-ri) 2 -30r 2 (u-r 2 ) 2 
-30r 2 (u-ri) 2 -30r^(u-r 2 ) 2 



A = J (i/ + z l ) dm , 



D 



(z 2 + x 2 ) dm , 



(123) 



+3r 4 + 3rf + 6r(r 2 



dmi dm 2 



(120) 

where all linear terms in ri or r 2 have been omitted since 
these two vectors are expressed relative to the barycen- 
ter of the respective body and their integral vanishes. In 
the section |2.3[ an additional hypothesis is made on the 
mass distribution of each body that simplifies the poten- 
tial. They are supposed to be symmetrical relative to the 
planes perpendicular to the principal axes of inertia. As 
a consequence, the integral of any odd power of ri or r 2 
cancels. 

Appendix B. Inertia integral 

Inertia integrals of homogeneous ellipsoids are computed 
in the following way. Let a, b, c be the three semi-axes of 
a homogeneous ellipsoid £ of density p. The total mass 
of the ellipsoid is 

4-7T 

to = —pabc (121) 
and the second order inertia integrals are 



px dxdydz = pa bc / X dX dY dZ = —ma , 
e Jb 5 

/ py 2 dxdydz = pab 3 c [ Y 2 dX dY dZ = -mb 2 , 
Js Jb 5 

f 1 

pz 2 dx dy dz = pabc 3 / Z 2 dX dY dZ = -mc 2 , 

s Jb 5 



(122) 



C = J (x 2 + y 2 )dm 



it: 

we get relations between the semi-axes and the moments 
of inertia 

a 2 = J^(-A + B + C) , 



2m 

b2 = i^ {A ~ B+c) 



(124) 



2to 



(A + B-C) 



General expressions of the inertia integrals are thus 

X i Y j Z k dXdYdZ 



x i y J z dm — - — ma l b 3 c k 



3 

47r' 



(125) 



with a, b, c given by (124 1 



Appendix C. Averaged quantities 

In this appendix, we give general formulas for the averag- 
ing over the orbital mean motions. The integrals will be 
computed using the true anomaly (y) as an intermediate 
variable. We recall first the basic formulas 



dM 
X 

y 



r 



a 2 ^/T— 
r cos v 
rsinv 
a{l~e 2 ) 
1 + e cos v 



-dv 



(126) 



where X and y are the coordinates of a point on a ke- 
plerian orbit in the reference frame (i,j,k) with i and 
k respectively in the direction of periapse and angular 
momentum. 

Intermediate integrals 

In the following, we handle integrals such as Wallis inte- 
grals. We recall their expression. Let 

-i /"27T -i p2tt 

I n =— cos 11 tdt=— sin n tdt (127) 
2vr J Q 2n J Q 
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and 



J 7i, rn — 



2tt 



27T 



sin™ ^ cos" i (ft 



2tt Jo 



2tt 



(128) 



sin" i cos™ i (ft. 



We have then 



l 2 2 P(p!) 2 



if n is odd, 

if n = 2p, p e N , 



(129) 



and for p > 0, ^2(p+i) can be computed using the recur- 
rence formula 

J 2(p+i) = ^t^ /2 p- ( 13 °) 

The integrals J„ iTO are null whenever n or m is odd, 
else their values are a sum of integrals Ik 

1 f 27T 

■h v 2q = 7z- / (1 - cos 2 t) p cos 29 i eft, 

= t(-D fc (; 

fe=0 v 

=X>1)*(0 



hq+2k 



(131) 



'2p+2fc- 



The last equality comes from J„ )TO = J m , n . 



Computation of (1/r") for n > 2 
From these results, we can write 

/■27T 

r^= ■ . '.. , : I 

, r 

E(n/2-l) 



1 \ 1 1 /" 27r 

— > = — ^ ^ — tt^ — / (1 + ecosv) n ~ 2 dv, 

r n I a n (l - e 2 )™- 3 / 2 2tt J y ' 



i ^ / n — 2\ 2p 

a«(l - e 2 )"-3/ 2 ^ I 2p ) i2p6 ' 

v 7 p— 



a"(l - e 2 )™- 3 / 2 



E(n/2-l) 



p=0 



where E(x) returns the integer part of x and 

/n-2\ (2p)l 
^" W ^ 2p ) 2 2 p( p !) 2 ' 

The recurrence relation for A n {p) is 



(132) 
(133) 



^(P+1)= ( "" 2P ( 2 P % 2P " 3) ^(P) ( 134 ) 
with A»(0) = 1. 



Computation of (X m y n /r l ) for Z > to + n + 2 

In averaging computations we meet integrals in the form 
' (r- si) kl ■ ■ ■ (r- Sj)^' 

These integrals can be computed from 
' x m y n \ _ 1 1 

a h(l- e 2)fc-3/2 2^ 



(135) 



/ ' 



x / cos™ i^sin" + ecosu) h 2 dv 



h-2 

I ^ /h-2 



a h (l - e 2 )^-3/2 V 2k 



fe=0 

(136) 

where h = l — m — n and J„ jTO defined as previously. This 
integral is null whenever n is odd. 
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